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We present an implementation of the adiabatic spin-wave 
dynamics of Niu and Kleinman. This techniques allows to 
decouple the spin and charge excitations of a many-electron 
system using a generalization of the adiabatic approximation. 
The only input for the spin-wave equations of motion are the 
energies and Berry curvatures of many-electron states describ- 
ing frozen spin spirals. The latter are computed using a newly 
developed technique based on constrained density-functional 
theory, within the local spin density approximation and the 
pseudo-potential plane-wave method. Calculations for iron 
show an excellent agreement with experiments. 



Over the last 30 years, Density-Functional Theory 
(DFT) ^ has brought a large number of properties of real 
materials within the predictive range of first-principles 
calculations. Even though excited states cannot be prop- 
erly described within DFT — which strictly speaking only 
applies to the electronic ground states — the use of the 
adiabatic approximation allows to predict the vibrational 
excitation spectra of semiconductors ||^ and metals ||] 
to a very high accuracy. Magnetic excitations are also 
adiabatically decoupled from charge excitations. How- 
ever, the adiabatic decoupling of spin excitations has 
only been implemented so far using ad-hoc models — such 
as the Heisenberg Hamiltonian ^ — which are difficult 
to reconcile with an itinerant picture of magnetism. A 
substantial step toward the calculation of magnetic ex- 
citations within a firm theoretical framework has been 
done by Niu and Kleinman (NK) who showed how 
a generalization of the Born-Oppenheimer (BO) method 
can be used to rigorously, although approximately, de- 
couple spin and charge excitations of a many-electron 
system. The NK method lends itself quite naturally to 
be implemented within DFT. The latter requires a proper 
description of non-coUinear magnetic structures, a task 
which has hardly been tackled so far in full generality. 
Although the local spin-density approximation (LSDA) 
was formulated long ago so as to account for general mag- 
netic structures the vast majority of the LSDA cal- 
culations assumes that the magnetization is aligned to a 
same direction all over the system. Most of the existing 
applications to non coUinear structures rely on some sort 
of atomic-sphere approximation (ASA) in which different 
s^pin quantization axes are chosen within different spheres 
|[7|. Progresses beyond this rather crude approximation 
have been made only recently ]^-[lO[| . 



In the present work we combine the adiabatic decou- 
pling method of Ref. |^ with a newly developed con- 
strained DFT scheme which allows the calculation of the 
spin-spiral states necessary to implement the NK spin- 
wave equations of motion, without using any atomic- 
sphere approximation. Our implementation — which is 
based on the plane-wave pseudo-potential method — is 
successfully demonstrated by calculating the magnon dis- 
persions of iron. This paper is organized as follows: we 
first briefly review and generalize the NK adiabatic de- 
coupling scheme; we then introduce our constrained DFT 
approach to spin-spiral states, and we finally present our 
results for iron. 

In order to establish an adiabatic decoupling scheme 
that allows to derive from first principles the spin- 
fiuctuation equations of motion, we first derive the clas- 
sical equations of motion for nuclei in a molecule or in a 
solid by using of the concept of adiabatic coherent state. 
This procedure has the advantage of not making explicit 
reference to the smallness of the nuclear kinetic energy 
term in the Hamiltonian. Because of this, we can then 
make the ansatz that a similar method applies in the gen- 
eral case where no such small term exists in the Hamilto- 
nian, but the density of states of the system is still char- 
acterized by a low-lying portion which is well separated 
in energy from the rest. This procedure provides thus an 
unified approach to the classical BO approximation and 
to the NK method. 

Let us consider a system of interacting electrons (of 
mass TO = 1) and nuclei (of mass M 3> to and in number 
A^), and let us indicate by capital letters the coordinates 
and momenta of the nuclei, {R,P), and by lower-case 
letters those of the electrons, {r,p). The Hamiltonian of 
the system reads: 
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p2 + V{r, R) 



(1) 



(here and in the following quantum operators are in- 
dicated by boldfaces). Let us now consider the wave- 
function ^'7^,•p (r, R) which minimizes the expectation 
value of the Hamiltonian, Eq. (^, with the constraint 
that the expectation values of R and P equal TZ and V, 
respectively. By introducing the A and fi Lagrange multi- 
pliers for the constraints on R and P respectively, "^iz.-p 
can be formally obtained as the ground-state solution of 
the eigenvalue equation: 

[H - X{R -TZ)- ^iiP - V)] ^n.v = P)^'7^,P, (2) 
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where A = A (7?., V) and /i = /x(7?., V) are chosen in such a 
way that: {'^n,v\R\'^n,v) = and {^-fi,v\P\'^n.v) = 
V. It is easy to verify that when H is the Hamiltonian of 
a harmonic osciUator, H — ^^P'^ + ^kR^, the eigenvalue 
equation, Eq. IM) yields directly Glauber's canonical 
coherent states [0. For this reason, in the general case 
we name ^'TC,-p(r, i?) an adiabatic coherent state (ACS) 
[pl and we refer to (7^, V) as to the label of the ACS. 

Let us suppose that the time evolution of an ACS is an 
ACS: this is true for harmonic oscillators and canonical 
coherent states, but in general this is only an approxima- 
tion. The low energy dynamical properties of the system 
are then given by the variational problem: 
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«'K(t),P(t) }dt = 0, (3) 



which, following the NK procedure, leads the equations 
of motion for the {TZ, V) label: 
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where the 
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ri's are 



matrices defined as Vif'-'^ — 



dAi 



and the i or j subscripts indicate 



the components of the TZ and V 3A^-dimensional vectors. 
From this definition, it follows that Vl- '- — —Vl- . Fol- 
lowing the literature , we name the Q. matrices the 
Berry curvature of the ACS. In order to recover the clas- 
sical equations of motion for the nuclear degrees of free- 
dom, we first simplify the eigenvalue equation (H) using 
the fact that H is quadratic in P. Taking into account 
that [P, e*'''-^ — V e^^^, a simple algebraic manipula- 
tion in Eq. (H) allows to cast '^^n.v and e^R^V) into the 
form: 



*7j,p(r,i?) =e'^^$7j(r,i?). 



(5) 
(6) 



where the real wave-function $7^ satisfies the eigenvalue 
equation: 



[^^-A(i^-7^)]$7^ = e(7^)$7^, 



(7) 



and the Lagrange multiplier A = \{TZ) is chosen in such 
a way that: {<^ti\R\^ti) = TZ. Using Eq. (||) and the re- 
ality of the il's can be explicitly calculated: fi,- 



en. 



= = and 

J ^ . -J ^ 0. By inserting these fi's and the expres- 
sion given by Eq. for e{TZ, V) into Eq. (^, one finally 
arrives at a set of Hamilton-like equations for the nuclear 
motion. In order identify e{TZ) with the familiar BO en- 
ergy surface, it is enough to neglect in Eq. (|^) the nuclear 
kinetic-energy term present in w 0). It is impor- 

tant to remark that — while the smallness of 1/Af is the 



physical reason why nuclear vibrations are adiabatically 
decoupled from electronic excitations — this property has 
not been formally used to establish the classical equations 
of motion, but only to identify e{TZ) with the traditional 
BO energy surface. In principle, e{TVi could be as well 
calculated from its definition, Eq. (|^), without letting 
1/M go to 0. When studying the electronic spectrum 
of magnetic systems, no such small term exists in the 
Hamiltonian. Nevertheless, spin excitations can be adia- 
batically decoupled from charge excitations by using the 
magnetic BO surface obtained from the magnetic analog 
ofEq. ^. 

Spin excitations are characterized by the fluctuation 
of the spin polarization, M{f), much in the same way as 
nuclear vibrations are characterized by the fluctuation 
of the nuclear canonical coordinates, {TZ^V). In both 
cases, the fluctuating observable can be identified with 
the label of an appropriate ACS. Magnetic ACS (MACS) 
are thus labeled by a vector field M{r)^ ^[MY ^'^^ ^^^■^ 
are solution of the variational problem: 



E[M] = min(*|J?|1'), 
with the constraint: 

(*|'0(f) + CT'0(f)|*) = 7W(f), 



(8) 



(9) 



where H is the many-body Hamiltonian of the sys- 
tem, the i/i's are fermion spinor field operators, and 
a = (da;, ay,(Tz) the Pauli matrices. 

Following NK, the linearized equations of motion for 
M{f,t) read {h = 1): 

^ / d\ "n^p{r\r")Mp{r",t) = 

J2 I d'r "K^f,{f\f")Mp{f",t), (10) 

where f2 is the Berry curvature of the MACS and K the 
matrix of the second derivatives of its energy with respect 
to its label, A4: 
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(12) 



In order to calculate the Berry curvature, we use a 
finite-difference approach ||,|l^. To keep the notation 
as simple as possible, we define ^'o = '^[M.ai'r ')], 
*i = ^[Ma{r) + i5ap5{r - r ')], and *2 = '^[Mo,{r) + 
eSa-y^ir - r ")]. We have then: VL^pir ',f ") w 
Im(*i-*o|^'2-*o) = Im((*i|5'2) + (*2|^'o) + 
(vJ^ol^i)). By using the relation: (*i|^'2) log(^'i 1*2) + 
1, we finally arrive at the relation: U,ap{r ',r ") = 
-^Im log ((«'o|4'i)(4'i|4'2)(*2|*o)). 

Before proceeding further, we transform the con- 
strained variational problem, Eqs. (j^,^), into a non- 
constrained one by making use of a penalty functional. 
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P[^], which favors those states whose magnetization, 
M[^]ir) = (^'|i/'(f)+CTi/)(r)|5'), is very close to the MACS 

label, M{r): P[*,A^] ^ A J {Mi^]{r) - M{r)fd^r. In 
the limit where A +00, the state which minimizes 
the functional: E[M,^ = + A /(M^.(f) - 

A4{r))^d'^r without any constraints but the trivial one 
on state normalization, will solve the variational prob- 
lem, Eq. (||), with the constraint given by Eq. Vari- 
ation with respect to leads to the eigenvalue equation: 

Jf* + 2AY^ct, ■ (A%](rO - Min))^ = 

i 

Ea[M\^, (13) 

where the sum runs over electrons. 

Eq. (^3|) is formally equivalent to a many-body 
Schrodinger equation where, in addition to physical in- 
teractions, electrons are subject to an external magnetic 

field whose magnitude, B{r) = 2yl ^M[^] (r) — A^(f)^ 
depends self-consistently upon the solution of the equa- 
tion, 4". When A grows large, B has the effect to 
make A/[,j,](r) closer and closer to A^(r), so that B{r) = 

2A(M[^](f) — M{r) tends to a finite hmit, and the value 
of the penalty functional tends to zero. We conclude 
that in the large- A limit the ground-state energy of the 
Schrodinger equation, Eq. (|l3), gives the value of the 
constrained minimum, Eqs. (p,^): Vm\A^+oo Ea[M\ = 
E[M\. 

An approximate solution of Eq. ( p^ can be found 
using DFT within the LSDA, provided that the latter 
is formulated in such a way as to allow for arbitrary, 
non coUinear, magnetic structures. In the LSDA, the 
many-body Hamiltonian appearing in Eq. ( p^ is re- 
placed with an independent-electron Hamiltonian where 
the effective one-body potential depends self-consistently 
upon the electron charge density and spin polarization. 
The implementation of the constraining penalty function 
in a non-coUinear LSDA code thus requires a trivial mod- 
ification of the dependence of the self-consistent potential 
upon the spin polarization. 

Inspection of the microscopic distribution of the mag- 
netization resulting from various LSDA calculations 
shows that it is non-vanishing mostly in the neighbor- 
hood of atoms (or molecules, in magnetic molecular crys- 
tals), and that within these neighborhoods its direc- 
tion remains almost constant even when different atomic 
(molecular) moments are not aligned. At variance with 
the assumptions which underlie the ASA treatment of 
(non-coUinear) magnetism, small deviations from perfect 
alignment may occur even within an atomic (or molec- 
ular) volume, and this quasi-coUinearity results from 
the contribution of many occupied Kohn-Sham orbitals 
whose individual magnetization is far from coUinear. In 
view of this fact, we choose to label the MACS with 
atomic (or molecular) moments, rather than with a con- 
tinuous magnetization. This discretization introduces a 



natural coarse graining in the treatment of the excita- 
tion spectrum, by cutting away all the fluctuations whose 
wavelength is shorter than the inter-atomic (-molecular) 
distance, and whose high-energies would be incompati- 
ble with the adiabatic decoupling scheme adopted here. 
Atomic magnetic moments are defined through the rela- 
tion: rh{R) — J dr u;(r—-R) 7W(r), where the f^'s indicate 
atomic positions, and wif) is a weight function which is 
equal to one in a region around the origin and zero else- 
where. In our applications, we choose these regions to 
be the largest possible non-overlapping spheres. Within 
this picture, long-wavelength magnetic excitations can 
be thought as time-dependent deviations of the localized 
moments from their ground-state configuration. Analo- 
gously to the current terminology adopted for vibrational 
excitations, a pattern of atomic moments, {m(i?)} is 
called a frozen magnon, and the corresponding matrix of 
second derivatives of the energy, d^E/dma{R)dm0{R'), 
the matrix of the inter- atomic torque constants. The self- 
consistent equations which yield the MACS correspond- 
ing to a frozen magnon have a form similar to that of Eq. 
(O) , where the self-consistent constraining magnetic field 
is different from zero only within the atomic volumes. For 
any given wave-vector of the magnon, this equation can 
be solved using super-cell techniques, much in the same 
way as these are applied to phonon calculations. 

If one restricts frozen magnons to configurations de- 
scribed by magnetic spirals | p5[ |, then the corresponding 
MACS can be calculated without any super-cell, by using 
a generalized Bloch theorem . The magnetic con- 

figuration of a spiral of wave- vector q is such that atomic 
moments of fixed magnitude form an angle 9 with respect 
to the z direction, independent of R, while the polar an- 
gle of rn{R) has the form: 0(^) = 0(0) +q- R. 

In order to work out the spin-wave equations of mo- 
tion, let us consider the special case of a mono-atomic 
ferro-magnetic crystal, such as bcc iron. In our ap- 
plications, we constrain 0, while the magnitude of the 
atomic moments is left unconstrained, and it turns out 
to be independent of 9 for small angles, thus giving an 
a-posteriori justification of the spiral model. Transla- 
tional invariance implies that Q. and K depend on R and 
R' only through their differences; rotational invariance 
around the z axis implies that K^y = Kyx = while 
the fact that = ^yy — follows from their defi- 
nition, Eq. (|Tl|). Using these relations, the equation 
of motion for the localized moments can be cast into 
the form: -iCl^yiq) ^"^(9) = Kxx{q)rh{q), where A{q) 

indicates the Fourier transform of A{R), and m{q) = 
mx{q) + irhyiq)- Magnon frequencies are thus given by: 

uj{ci) = kxx{q)/^xy{q)- 

We have applied this method to the calculation of 
magnon frequencies in bcc iron. The computations are 
performed using the LSDAjT^ within the plane-wave 
pseudo-potential method [|18[. Plane waves up to a 
kinetic-energy cutoff of 90 Ry have been included in 
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the basis set, while Brillouin-zone integrations were per- 
formed on a 8 X 8 X 8 mesh, with a Gaussian-smearing 
of 0.03 Ry In Fig. 1 we display the results of our 

calculations, performed for q along the (100) direction, 
while the ground-state magnetization axis was chosen to 
be (001). The agreement with available experiments is 
very good, thus giving confidence in the general theoret- 
ical framework followed in this work. 

We believe that the adiabatic decoupling scheme by 
Niu and Kleinman, based on the concept of adiabatic co- 
herent state, still needs stronger theoretical foundations. 
However, we have demonstrated that this scheme nat- 
urally leads to the well established Born-Oppenheimer 
classical equation of motion for nuclei in molecules or 
solids. Furthermore, this scheme can be straightfor- 
wardly implemented within density-functional theory 
and, when applied to magnetic excitations of bcc iron, 
it leads to magnon dispersions in excellent agreement 
with experiment. We conclude that further theoretical 
and computational work will probably demonstrate the 
effectiveness of this method to deal with slow dynamics 
in quantum systems. 
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FIG. 1. Calculated and experimental magnon dispersions 
in iron. Open hexagons: pure iron at lOK [20]. Open tri- 
angles: Fe with 12 % Si) at room temperature [21]. Full 
hexagons: this work. The line is a guide for the eye. 
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